#include "aerosol_atmosphere.h"
#include <QDebug>
#include <QFile>
#include <QString>
#include <QDir>
#include <QDateTime>
#include "mie_c_src/mie.h"

double (*AerosolAtmosphere::m_LSameDir)[N_ZENITH][2*N_PI][2][2];
double (*AerosolAtmosphere::m_LDiffDir)[N_ZENITH][2*N_PI][2][2];

int (*AerosolAtmosphere::m_scatAnglesIndexSameDir)[N_ZENITH][2*N_PI];
int (*AerosolAtmosphere::m_scatAnglesIndexDiffDir)[N_ZENITH][2*N_PI];

AerosolAtmosphere::AerosolAtmosphere(QString layerName,
				     double hBottom,
				     double hTop,
				     double (*smatrixAero)[3][3],
				     double (*pbsdfAeroAveraged)[N_ZENITH][2*N_PI][3][3],
				     double aerosolTau,
				     QObject *parent) :
    AtmosLayer(layerName, hBottom, hTop, parent),
    m_aeroTau(aerosolTau),
    m_smatrixAero(smatrixAero),
    m_averagedFowardPbsdf(pbsdfAeroAveraged)
{
    if(m_smatrixAero == 0){
	qDebug() << "No valid buffer for aerosol scattering matrix provided.";
	exit(0);
    }

    m_smatrixAeroAndMole = new double[N_SCAT][3][3];
    for(int i = 0; i < N_SCAT; i++){
	for(int j = 0; j < 3; j++){
	    for(int k = 0; k < 3; k++){
		m_smatrixAeroAndMole[i][j][k] = 0;
	    }
	}
    }

    m_tau = m_moleTau + m_aeroTau;
    scatMatrixCaculate();

    if(m_LSameDir == 0){
	LMatrixInit();
    }
    if(m_scatAnglesIndexSameDir == 0){
	scatAngleIndexInit();
    }
}

AerosolAtmosphere::~AerosolAtmosphere()
{
    delete []m_smatrixAeroAndMole;
    if(m_scatAnglesIndexSameDir != 0){
	delete []m_scatAnglesIndexSameDir;
	delete []m_scatAnglesIndexDiffDir;
	delete []m_LSameDir;
	delete []m_LDiffDir;
    }
}

//#define LMATRIX_DEBUG
void AerosolAtmosphere::LMatrixInit()
{
    m_LSameDir = new double[N_ZENITH][N_ZENITH][2*N_PI][2][2];
    m_LDiffDir = new double[N_ZENITH][N_ZENITH][2*N_PI][2][2];

//    double muScat = 0;
    double thetaI, phiI, thetaS, phiS;
    thetaI = phiI = thetaS = phiS = 0;

    //same zenith direction
    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	thetaI = m_zenithUpList[iThetaI];
	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    thetaS = m_zenithUpList[iThetaS];
	    for(int iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

		if(iThetaI == iThetaS && iPhiS == 0){
		    //complete same dir
		    m_LSameDir[iThetaI][iThetaS][iPhiS][0][0] = -1;
		    m_LSameDir[iThetaI][iThetaS][iPhiS][0][1]= 0;
		    m_LSameDir[iThetaI][iThetaS][iPhiS][1][0]= -1;
		    m_LSameDir[iThetaI][iThetaS][iPhiS][1][1]= 0;
		    continue;
		}

		phiS = m_aziList[iPhiS];

		m_LSameDir[iThetaI][iThetaS][iPhiS][0][0] = (-(Power(Sin(phiI - phiS),2)*Power(Sin(thetaS),2)) + Power(Cos(thetaS)*Sin(thetaI) - Cos(phiI - phiS)*Cos(thetaI)*Sin(thetaS),2))/(1 - Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2));
		m_LSameDir[iThetaI][iThetaS][iPhiS][0][1] = (2*Sin(phiI - phiS)*Sin(thetaS)*(-(Cos(thetaS)*Sin(thetaI)) + Cos(phiI - phiS)*Cos(thetaI)*Sin(thetaS)))/(1 - Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2));
		m_LSameDir[iThetaI][iThetaS][iPhiS][1][0] = (-(Power(Sin(phiI - phiS),2)*Power(Sin(thetaI),2)) + Power(Cos(phiI - phiS)*Cos(thetaS)*Sin(thetaI) - Cos(thetaI)*Sin(thetaS),2))/(1 - Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2));
		m_LSameDir[iThetaI][iThetaS][iPhiS][1][1] = (2*Sin(phiI - phiS)*Sin(thetaI)*(Cos(phiI - phiS)*Cos(thetaS)*Sin(thetaI) - Cos(thetaI)*Sin(thetaS)))/(1 - Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2));

#ifdef LMATRIX_DEBUG
		if(iThetaI == 6 && iThetaS == 20 && iPhiS == 40){
		    qDebug() << thetaI << thetaS << phiS
			     << m_LSameDir[iThetaI][iThetaS][iPhiS][0][0]
			     << m_LSameDir[iThetaI][iThetaS][iPhiS][0][1]
			     << m_LSameDir[iThetaI][iThetaS][iPhiS][1][0]
			     << m_LSameDir[iThetaI][iThetaS][iPhiS][1][1];
		}
#endif
	    }
	}
    }

    //different zenith direction
    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	thetaI = m_zenithUpList[iThetaI];
	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    thetaS = m_zenithDList[iThetaS]; //down, different
	    for(int iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

		if(iThetaI == iThetaS && iPhiS == N_PI){
		    //complete reverse dir
		     m_LDiffDir[iThetaI][iThetaS][iPhiS][0][0] = -1;
		     m_LDiffDir[iThetaI][iThetaS][iPhiS][0][1]= 0;
		     m_LDiffDir[iThetaI][iThetaS][iPhiS][1][0]= -1;
		     m_LDiffDir[iThetaI][iThetaS][iPhiS][1][1]= 0;
		    continue;
		}

		phiS = m_aziList[iPhiS];

		m_LDiffDir[iThetaI][iThetaS][iPhiS][0][0] = (-(Power(Sin(phiI - phiS),2)*Power(Sin(thetaS),2)) + Power(Cos(thetaS)*Sin(thetaI) - Cos(phiI - phiS)*Cos(thetaI)*Sin(thetaS),2))/(1 - Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2));
		m_LDiffDir[iThetaI][iThetaS][iPhiS][0][1] = (2*Sin(phiI - phiS)*Sin(thetaS)*(-(Cos(thetaS)*Sin(thetaI)) + Cos(phiI - phiS)*Cos(thetaI)*Sin(thetaS)))/(1 - Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2));
		m_LDiffDir[iThetaI][iThetaS][iPhiS][1][0] = (-(Power(Sin(phiI - phiS),2)*Power(Sin(thetaI),2)) + Power(Cos(phiI - phiS)*Cos(thetaS)*Sin(thetaI) - Cos(thetaI)*Sin(thetaS),2))/(1 - Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2));
		m_LDiffDir[iThetaI][iThetaS][iPhiS][1][1] = (2*Sin(phiI - phiS)*Sin(thetaI)*(Cos(phiI - phiS)*Cos(thetaS)*Sin(thetaI) - Cos(thetaI)*Sin(thetaS)))/(1 - Power(Cos(thetaI)*Cos(thetaS) + Cos(phiI - phiS)*Sin(thetaI)*Sin(thetaS),2));

#ifdef LMATRIX_DEBUG
		if(iThetaI == 6 && iThetaS == 20 && iPhiS == 40){
		    qDebug() << thetaI << thetaS << phiS
			     << m_LDiffDir[iThetaI][iThetaS][iPhiS][0][0]
			     << m_LDiffDir[iThetaI][iThetaS][iPhiS][0][1]
			     << m_LDiffDir[iThetaI][iThetaS][iPhiS][1][0]
			     << m_LDiffDir[iThetaI][iThetaS][iPhiS][1][1];
		}
#endif

	    }
	}
    }
}

//#define MASK_FOWARD
void AerosolAtmosphere::scatMatrixCaculate()
{
    double moleRatio, aeroRatio;
    moleRatio = m_moleTau/(m_moleTau + m_aeroTau);
    aeroRatio = m_aeroTau/(m_moleTau + m_aeroTau);
    m_moleRatio = moleRatio;
    m_aeroRatio = aeroRatio;

    double angleList[N_SCAT];
    double angle = 0;
    double mu = 1;
    for(int i = 0; i < N_SCAT; i++){
	mu = cos(angle);
	angleList[i] = angle;

	m_smatrixAeroAndMole[i][0][0] = aeroRatio * m_smatrixAero[i][0][0] + 0.75*(1 + mu*mu)*moleRatio;
	m_smatrixAeroAndMole[i][0][1] = aeroRatio * m_smatrixAero[i][0][1] + 0.75*(-1 + mu*mu)*moleRatio;
	m_smatrixAeroAndMole[i][1][0] = aeroRatio * m_smatrixAero[i][1][0] + 0.75*(-1 + mu*mu)*moleRatio;
	m_smatrixAeroAndMole[i][1][1] = m_smatrixAeroAndMole[i][0][0];
	m_smatrixAeroAndMole[i][2][2] = aeroRatio * m_smatrixAero[i][2][2] + 1.5*mu*moleRatio;

	if(i < N_SMALL_SCAT){
	    angle += STEP_SMALL_SCAT;
	} else{
	    angle += STEP_BIG_SCAT;
	}
    }

}

//#define SCAT_ANGLE_INDEX_DEBUG
void AerosolAtmosphere::scatAngleIndexInit()
{
#ifdef SCAT_ANGLE_INDEX_DEBUG
    qDebug() << "AerosolAtmosphere::scatAngleIndexInit()";
#endif
    if(m_zenithUpList.isEmpty()){
	zenithAngleInit();
	solidAngleCaculate();
	scatAngleInit();
    }

    m_scatAnglesIndexSameDir = new int[N_ZENITH][N_ZENITH][2*N_PI];
    m_scatAnglesIndexDiffDir = new int[N_ZENITH][N_ZENITH][2*N_PI];

    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    for(int iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){
		int index = 0;
		double scatAngle = 0;
		//same zenith direct
		scatAngle = m_scatAnglesSameDir[iThetaI][iThetaS][iPhiS];
		if(scatAngle < SMALL_SCAT_REGION){
		    index = qRound(scatAngle / STEP_SMALL_SCAT);
		} else{
		    index = qRound((scatAngle-SMALL_SCAT_REGION) / STEP_BIG_SCAT + N_SMALL_SCAT - 1);
		}
		m_scatAnglesIndexSameDir[iThetaI][iThetaS][iPhiS] = index;
#ifdef SCAT_ANGLE_INDEX_DEBUG
		if(iThetaI == 7 && iThetaS == 20 && iPhiS == 40)
		    qDebug() << iThetaI << iThetaS << iPhiS << m_scatAnglesSameDir[iThetaI][iThetaS][iPhiS] << m_scatAnglesIndexSameDir[iThetaI][iThetaS][iPhiS];
#endif
		//different zenith direct
		scatAngle = m_scatAnglesDiffDir[iThetaI][iThetaS][iPhiS];
		if(scatAngle < SMALL_SCAT_REGION){
		    index = qRound(scatAngle / STEP_SMALL_SCAT);
		} else{
		    index = qRound((scatAngle-SMALL_SCAT_REGION) / STEP_BIG_SCAT + N_SMALL_SCAT - 1);
		}
		m_scatAnglesIndexDiffDir[iThetaI][iThetaS][iPhiS] = index;
#ifdef SCAT_ANGLE_INDEX_DEBUG
		if(iThetaI == 7 && iThetaS == 20 && iPhiS == 40)
		    qDebug() << iThetaI << iThetaS << iPhiS << m_scatAnglesDiffDir[iThetaI][iThetaS][iPhiS] << m_scatAnglesIndexDiffDir[iThetaI][iThetaS][iPhiS];
#endif
	    }
	}
    }

}

//#define DEBUG_LogAeroCaculate
void logAeroCaculate(double wl, double *ratio,
		     double rMin, int nRStep,
		     struct complex m1,
		     struct complex m2,
		     double *mu,
		     double pbsdf[][3][3], double *crossAvg);

void logDistribution(double rMean, double sigma,
		     double *ratio, double rMin, long nStep);
/*! \~english
    \fn AbstractSerial::AbstractSerial(QObject *parent)
    Averaged aerosol scattering matrix at N_SCAT angles for a two peak Log-Normal distribution.\n
    \note The subscript 1 is for small particals, 2 for big ones. And, you must generate a buffer m_smatrixAero[N_SCAT][3][3] before call this function, and pass the pointer to it, output will save in buffer for latter use.
*/
void AerosolAtmosphere::aerosolInitTwoPeak(double rMean1,
			double sigma1,
			struct complex m1,
			double rMean2,
			double sigma2,
			struct complex m2,
			double numRatio1_2, //tau ratio for small/big
			double (*smatrix)[3][3])
{
    if(rMean1 > R_MAX || rMean1 < R_MIN ||
	    rMean2 > R_MAX || rMean2 < R_MIN){
	    qDebug() << "Aerosol partical radius must between 0.01 and 10, while you input rMean:" << rMean1 << rMean2;
	    exit(82); //ASCII code for R
    }

    QString fileNamePrefix("scatMatrix/");
//    fileNamePrefix.append(QDateTime::currentDateTime().toString("MMDD_HH_MM_SS"));
    fileNamePrefix.append(QString().sprintf("R_%.2f_%.2f-M1_%.4f_%.4f-M2_%.4f_%.4f-Ratio_%.2f_Width_%.2f_%.2f", rMean1, rMean2, m1.re, m1.im, m2.re, m2.im, numRatio1_2, sigma1, sigma2));
    QDir dir(QDir::current());
    dir.mkpath(fileNamePrefix);
    if(QFile(fileNamePrefix+"/aerosol_scat_matrix.txt").exists()){
        qDebug() << "AerosolAtmosphere::aerosolInitTwoPeak(): current aerosol's micro-physical state has a scattering matrix template:" <<
                    QString(fileNamePrefix+"/aerosol_scat_matrix.txt");
        readAerosolScatMatrix(QString(fileNamePrefix+"/aerosol_scat_matrix.txt"),
                              smatrix);
        return;
    }
/***************Caculate number desity profile****************/
    double *ratioSmall = new double[NUM_R_STEPS];
    double *ratioBig = new double[NUM_R_STEPS];
    double *ratio = new double[NUM_R_STEPS];
    for(int i = 0; i < NUM_R_STEPS; i++){
	ratioSmall[i] = ratioBig[i] = ratio[i] = 0;
    }
    logDistribution(rMean1, sigma1, ratioSmall, R_MIN, NUM_R_STEPS);
    logDistribution(rMean2, sigma2, ratioBig, R_MIN, NUM_R_STEPS);
    for(int i = 0; i < NUM_R_STEPS; i++){
	double part2 = 1/(1+numRatio1_2);
	double part1 = 1- part2;
	ratio[i] = part1 * ratioSmall[i] + part2 * ratioBig[i];
    }
    {//save aerosol size spectrum
	QByteArray data("Radius\tRatio\n");
	QString line;
	double r = R_MIN;
	for(int i = 0; i < NUM_R_STEPS; i++, r+=R_STEP){
	    line.sprintf("%f\t%f\t%f\t%f\n", r, ratio[i],ratioSmall[i],ratioBig[i]);
	    data.append(line);
	}
	qDebug() << "Saving size spectrum...";
	QFile file(fileNamePrefix+"/size_spectrum.txt");
	file.open(QIODevice::ReadWrite);
	file.write(data);
	file.close();
    }

    if(smatrix == 0){
	qDebug() << "No valid buffer for aerosol scattering matrix provided.";
    }

    double mu[N_SCAT];
    double angleList[N_SCAT];
    double angle = 0;
    for(int i = 0; i < N_SCAT; i++){
	mu[i] = cos(angle);
	angleList[i] = angle;
	if(i < N_SMALL_SCAT){
	    angle += STEP_SMALL_SCAT;
	} else{
	    angle += STEP_BIG_SCAT;
	}
    }
    double crossAvg = 0;
    /***** start caculate small/big particals' scatMatrix ******/
    logAeroCaculate(m_wl, ratio, R_MIN, NUM_R_STEPS,
		    m1, m2, mu,
		    smatrix, &crossAvg);

#ifdef MASK_FOWARD
    double (*pbsdf)[3][3] = new double[N_SCAT][3][3];

    double m11, m12, m22, m33;
    m11 = m12 = m22 = m33 = 0;

    for(int i = 0; angleList[i] < ANGLE_STEP/2; i++){
	//first caculate average value across [angle, angle + STEP/2]
	int j = 0;
	m11 = m12 = m22 = m33 = 0;
	for(j = i; angleList[j] < angleList[i] + ANGLE_STEP/2; j++){
	    m11 += (smatrix[j][0][0] + smatrix[j+1][0][0]) / 2;
	    m12 += (smatrix[j][0][1] + smatrix[j+1][0][1]) / 2;
	    m22 += (smatrix[j][1][1] + smatrix[j+1][1][1]) / 2;
	    m33 += (smatrix[j][2][2] + smatrix[j+1][2][2]) / 2;
	}
	double tmp;
	tmp = (angleList[i] + ANGLE_STEP/2 - angleList[j-1]) / (angleList[j] - angleList[j-1]);

	pbsdf[i][0][0] = (smatrix[j][0][0]*tmp + smatrix[j-1][0][0]*(1-tmp)) * (angleList[i] + ANGLE_STEP/2 - angleList[j-1]) / (ANGLE_STEP/2);
	pbsdf[i][0][1] = (smatrix[j][0][1]*tmp + smatrix[j-1][0][1]*(1-tmp)) * (angleList[i] + ANGLE_STEP/2 - angleList[j-1]) / (ANGLE_STEP/2);
	pbsdf[i][1][1] = (smatrix[j][1][1]*tmp + smatrix[j-1][1][1]*(1-tmp)) * (angleList[i] + ANGLE_STEP/2 - angleList[j-1]) / (ANGLE_STEP/2);
	pbsdf[i][2][2] = (smatrix[j][2][2]*tmp + smatrix[j-1][2][2]*(1-tmp)) * (angleList[i] + ANGLE_STEP/2 - angleList[j-1]) / (ANGLE_STEP/2);

	pbsdf[i][0][0] += m11 * STEP_SMALL_SCAT / (ANGLE_STEP/2);
	pbsdf[i][0][1] += m12 * STEP_SMALL_SCAT / (ANGLE_STEP/2);
	pbsdf[i][1][1] += m22 * STEP_SMALL_SCAT / (ANGLE_STEP/2);
	pbsdf[i][2][2] += m33 * STEP_SMALL_SCAT / (ANGLE_STEP/2);
    }
    //then set matrix to average value at small angles
    for(int i = 0; angleList[i] < ANGLE_STEP/2; i++){
	smatrix[i][0][0] = pbsdf[i][0][0];
	smatrix[i][0][1] = pbsdf[i][0][1];
	smatrix[i][1][0] = smatrix[i][0][1];
	smatrix[i][1][1] = pbsdf[i][1][1];
	smatrix[i][2][2] = pbsdf[i][2][2];

    }
#endif

    saveScatMatrix(fileNamePrefix+"/aerosol_scat_matrix.txt", smatrix);
    qDebug() << "AerosolAtmosphere::aerosolInitTwoPeak(): template saved in: " + fileNamePrefix ;

    //delete when caculate over.
    delete []ratio;
    delete []ratioSmall;
    delete []ratioBig;
}


//caculate aerosol scatter property according to giving size spectrum
void logAeroCaculate(double wl, double *ratio,
		     double rMin, int nRStep,
		     struct complex mSmall,
		     struct complex mBig,
		     double *mu,
		     double pbsdf[][3][3], double *crossAvg)
{
	//***********caculate averaged aerosol properties.
	double sigmaRatio = 0; //cross section at each radius mutiply this radius's log-Norm distribution ratio, intermediate usage.
	struct complex s1[N_SCAT], s2[N_SCAT];
	double qext, qsca, qback, g;

	double radius = rMin;
	double x = 0; //size parameter: x = 2 PI radius / wavelength;
	double parallel;
	double perpen;
	for(int i = 0; i < nRStep; i++, radius += R_STEP){
	    if(ratio[i] == 0)
		continue;

	    x = 2 * M_PI * radius / wl;

	    //caculate scatter matrix of all angle at current radius
	    if(radius < R_REFREC_INDEX_SEPERATE)
		Mie(x, mSmall, mu, N_SCAT,
		    s1, s2, &qext, &qsca, &qback, &g);
	    else
		Mie(x, mBig, mu, N_SCAT,
		    s1, s2, &qext, &qsca, &qback, &g);

	    sigmaRatio = M_PI * pow(radius, 2) * qext * ratio[i];
	    *crossAvg += sigmaRatio;
	    //qDebug() << radius << *crossAvg << sigmaRatio;

	    for(int j = 0; j < N_SCAT; j++){

		parallel = 4 * cnorm (s2[j]) / (x * x * qext);
		perpen = 4 * cnorm (s1[j]) / (x * x * qext);
/* Bohren, P65, (3.16) and miemain.c:
   phasefn[i] = (parallel[i] + perpen[i]) / 2;
*/
		pbsdf[j][0][0] += sigmaRatio *
			(parallel + perpen) / 2;
		pbsdf[j][0][1] += sigmaRatio *
			(parallel - perpen) / 2;
		//fixed 4-23, 2012
		pbsdf[j][2][2] += sigmaRatio *
			(4*(s1[j].re * s2[j].re +
			    s1[j].im * s2[j].im)/(x*x*qext));
	    }
	}

	for(int j = 0; j < N_SCAT; j++){
		pbsdf[j][0][0] /= *crossAvg;
		pbsdf[j][0][1] /= *crossAvg;
		pbsdf[j][1][0] = pbsdf[j][0][1];
		pbsdf[j][1][1] = pbsdf[j][0][0]; //fixed 2012-4-8, see(C.F.Bohren equation 4.77)
		pbsdf[j][2][2] /= *crossAvg;
	}
}

//Log-Normal distribution ratio caculate
void logDistribution(double rMean, double sigma,
		     double *ratio, double rMin, long nStep)
{
	double r = rMin;
	double whole = 0;
	for(int i = 0; i < nStep; i++, r += R_STEP){
		ratio[i] = exp(-1 * pow((log(r) - log(rMean)), 2) /
			       (2 * sigma*sigma)) / r*sigma;
		whole += ratio[i];
//		qDebug() << r << ratio[i];

	}
	//	qDebug() << whole;

	for(int i = 0; i < nStep; i++){
		ratio[i] = ratio[i] / whole;
//		qDebug() << ratio[i];
	}
//	qDebug() << "****************************";
}


#define SAME_DIR true
#define DIFF_DIR false
void AerosolAtmosphere::scatCoeff(int iThetaS, int iPhiS,
				  int iThetaI, int iPhiI,
				  double lightIncident[3],
				  linearPolar *lightOut,
				  int scatMode)
{
    bool isSameDir = (scatMode==UP_IN_UP_SCAT || scatMode==D_IN_D_SCAT) ? SAME_DIR : DIFF_DIR;
    int (*scatAnglesIndex)[N_ZENITH][2*N_PI] = (isSameDir) ? m_scatAnglesIndexSameDir: m_scatAnglesIndexDiffDir;

    double (*LMatrix)[N_ZENITH][2*N_PI][2][2] = (isSameDir) ? m_LSameDir: m_LDiffDir;

    iThetaI = qAbs(iThetaI);
    iThetaS = qAbs(iThetaS);

    double M00, M01, M22;
    double L00, L01, L10, L11;

    int phiIndex = iPhiS - iPhiI;
    if(phiIndex < 0){
	phiIndex = 2*N_PI + phiIndex;
    }

    if(isSameDir && m_scatAnglesSameDir[iThetaI][iThetaS][phiIndex] < F_AVERAGE_REGION){
	double mu = Cos(m_scatAnglesSameDir[iThetaI][iThetaS][phiIndex]);
	M00 = m_aeroRatio*m_averagedFowardPbsdf[iThetaI][iThetaS][phiIndex][0][0] + 0.75*(1 + mu*mu)*m_moleRatio;
	M01 = m_aeroRatio*m_averagedFowardPbsdf[iThetaI][iThetaS][phiIndex][0][1] + 0.75*(-1 + mu*mu)*m_moleRatio;
	M22 = m_averagedFowardPbsdf[iThetaI][iThetaS][phiIndex][2][2] + 1.5*mu*m_moleRatio;
    } else{

	int iScatAngle = scatAnglesIndex[iThetaI][iThetaS][phiIndex];
	M00 = m_smatrixAeroAndMole[iScatAngle][0][0];
	M01 = m_smatrixAeroAndMole[iScatAngle][0][1];
	M22 = m_smatrixAeroAndMole[iScatAngle][2][2];
    }
    L00= LMatrix[iThetaI][iThetaS][phiIndex][0][0];
    L01= LMatrix[iThetaI][iThetaS][phiIndex][0][1];
    L10= LMatrix[iThetaI][iThetaS][phiIndex][1][0];
    L11= LMatrix[iThetaI][iThetaS][phiIndex][1][1];

    lightOut->i = lightIncident[I]*M00 + L00*lightIncident[Q]*M01 + L01*lightIncident[U]*M01;
    lightOut->q = L10*lightIncident[I]*M01 + lightIncident[U]*(L01*L10*M00 + L00*L11*M22) + lightIncident[Q]*(L00*L10*M00 - L01*L11*M22);
    lightOut->u = -(L11*lightIncident[I]*M01) + lightIncident[U]*(-(L01*L11*M00) + L00*L10*M22) + lightIncident[Q]*(-(L00*L11*M00) - L01*L10*M22);

    lightOut->i *= (1/(4*M_PI));
    lightOut->q *= (1/(4*M_PI));
    lightOut->u *= (1/(4*M_PI));

//    if(iThetaI==0 && iPhiI==90 && iThetaS==0 && iPhiS==0)
//	qDebug() << iThetaS << iPhiS << iThetaI << iPhiI
//		 << "incident" << lightIncident[0] << lightIncident[1] << lightIncident[2]
//		 << "out" << lightOut->i << lightOut->q << lightOut->u
//		 << "M" << M00 << M01 << M22
//		 << "L"<< L00 << L01 << L10 << L11;
}

void AerosolAtmosphere::solarScatCoeff(int iThetaS, int iPhiS,
				       linearPolar *lightOut,
				       int scatMode)
{
    //incident light direction is solar flux direction.
    int iThetaI = qRound(-1.0 *
			 (M_PI-m_solarZen-ANGLE_STEP/2)/ANGLE_STEP);
    int iPhiI = qRound(m_solarAzi/ANGLE_STEP);
//    qDebug() << "AerosolAtmosphere::solarScatCoeff: " <<  m_zenithDList[-iThetaI]*180/M_PI;
    //incident light parameter is solar flux.
    double lightIncident[3];
    lightIncident[0] = m_solarIntensity ;
    lightIncident[1] = 0;
    lightIncident[2] = 0;

/* below: same code as in scatCoeff()*/
    bool isSameDir = (scatMode==UP_IN_UP_SCAT || scatMode==D_IN_D_SCAT) ? SAME_DIR : DIFF_DIR;
    int (*scatAnglesIndex)[N_ZENITH][2*N_PI] = (isSameDir) ? m_scatAnglesIndexSameDir: m_scatAnglesIndexDiffDir;

    double (*LMatrix)[N_ZENITH][2*N_PI][2][2] = (isSameDir) ? m_LSameDir: m_LDiffDir;

    iThetaI = qAbs(iThetaI);
    iThetaS = qAbs(iThetaS);

    double M00, M01, M22;
    double L00, L01, L10, L11;

    int phiIndex = iPhiS - iPhiI;
    if(phiIndex < 0){
	phiIndex = 2*N_PI + phiIndex;
    }
    int iScatAngle = scatAnglesIndex[iThetaI][iThetaS][phiIndex];

    if(isSameDir && m_scatAnglesSameDir[iThetaI][iThetaS][phiIndex] < F_AVERAGE_REGION){
	double mu = Cos(m_scatAnglesSameDir[iThetaI][iThetaS][phiIndex]);
	M00 = m_aeroRatio*m_averagedFowardPbsdf[iThetaI][iThetaS][phiIndex][0][0] + 0.75*(1 + mu*mu)*m_moleRatio;
	M01 = m_aeroRatio*m_averagedFowardPbsdf[iThetaI][iThetaS][phiIndex][0][1] + 0.75*(-1 + mu*mu)*m_moleRatio;
	M22 = m_averagedFowardPbsdf[iThetaI][iThetaS][phiIndex][2][2] + 1.5*mu*m_moleRatio;
    } else{
	M00 = m_smatrixAeroAndMole[iScatAngle][0][0];
	M01 = m_smatrixAeroAndMole[iScatAngle][0][1];
	M22 = m_smatrixAeroAndMole[iScatAngle][2][2];
    }
    L00= LMatrix[iThetaI][iThetaS][phiIndex][0][0];
    L01= LMatrix[iThetaI][iThetaS][phiIndex][0][1];
    L10= LMatrix[iThetaI][iThetaS][phiIndex][1][0];
    L11= LMatrix[iThetaI][iThetaS][phiIndex][1][1];

    lightOut->i = lightIncident[I]*M00 + L00*lightIncident[Q]*M01 + L01*lightIncident[U]*M01;
    lightOut->q = L10*lightIncident[I]*M01 + lightIncident[U]*(L01*L10*M00 + L00*L11*M22) + lightIncident[Q]*(L00*L10*M00 - L01*L11*M22);
    lightOut->u = -(L11*lightIncident[I]*M01) + lightIncident[U]*(-(L01*L11*M00) + L00*L10*M22) + lightIncident[Q]*(-(L00*L11*M00) - L01*L10*M22);
/* above: same code as in scatCoeff()*/

//    lightOut->i *= (1/(4*M_PI));
//    lightOut->q *= (1/(4*M_PI));
//    lightOut->u *= (1/(4*M_PI));
    lightOut->i /= 4;
    lightOut->q /= 4;
    lightOut->u /= 4;

//    qDebug() << iThetaS << iPhiS << iThetaI << iPhiI
//	     << "incident" << lightIncident[0] << lightIncident[1] << lightIncident[2]
//	     << "out" << lightOut->i << lightOut->q << lightOut->u
//	     << "M" << M00 << M01 << M22
//	     << "L"<< L00 << L01 << L10 << L11
//	     << iScatAngle;

}

void saveScatMatrix(QString fileName, double (*pbsdfAero)[3][3])
{
    QFile file(fileName);
    QByteArray fileData;

    fileData.append("#scatAngle\tM11\tM12\tM13\tM21\tM22\tM23\tM31\tM32\tM33\n");
    for(int i = 0; i < N_SCAT; i++){
	double angle = 0;
	if(i < N_SMALL_SCAT){
	    angle = i * STEP_SMALL_SCAT;
	} else{
	    angle = SMALL_SCAT_REGION + (i-N_SMALL_SCAT+1) * STEP_BIG_SCAT;
	}

	fileData.append(QString().sprintf("%f\t", angle));
	for(int j = 0; j < 3; j++){
	    for(int k = 0; k < 3; k++){
		fileData.append(QString().sprintf("%f\t", pbsdfAero[i][j][k]));
	    }
	}
	fileData.append("\n");
    }

    file.open(QIODevice::ReadWrite);
    file.write(fileData);
    file.close();
}

#include <QTextStream>
//tested ok, 4-10, 2012, 9:55
void AerosolAtmosphere::readAerosolScatMatrix(QString fileName, double (*pbsdfAero)[3][3])
{
    QFile file(fileName);
    if(!file.exists()){
	qDebug() << "Error: AerosolAtmosphere::readAerosolScatMatrix(), file " + fileName + " not exists!";
	exit(-1);
    }
    file.open(QIODevice::ReadOnly);
    QByteArray data(file.readAll());
    file.close();

    QTextStream strm(data);

    strm.readLine(); //pass first line
    QString str, line;
    QTextStream lineStrm;
    for(int i = 0; i < N_SCAT; i++){
	line.clear();
	line = strm.readLine();
	lineStrm.setString(&line, QIODevice::ReadOnly);
	str.clear();
	lineStrm >> str; //pass first angle word.
	for(int j = 0; j < 3; j++){
	    for(int k = 0; k < 3; k++){
		str.clear();
		lineStrm >> str;
		pbsdfAero[i][j][k] = str.toDouble();
//		qDebug() << pbsdfAero[i][j][k];
	    }
	}
    }
}


#define STEP_AVG ANGLE_STEP/40
void AerosolAtmosphere::averageAerosolMatrix(double (*smatrix)[3][3], double (*pbsdfAero)[N_ZENITH][2*N_PI][3][3])
{
    if(m_scatAnglesIndexSameDir == 0){
	scatAngleIndexInit();
    }

    double phiI = 0;
    for(int iThetaI = 0; iThetaI < N_ZENITH; iThetaI++){
	double thetaI;
	thetaI = m_zenithUpList[iThetaI];

	for(int iThetaS = 0; iThetaS < N_ZENITH; iThetaS++){
	    double thetaS;
	    thetaS = m_zenithUpList[iThetaS];

	    for(int iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){
		double phiS;
		phiS = m_aziList[iPhiS];

		if(m_scatAnglesSameDir[iThetaI][iThetaS][iPhiS] < F_AVERAGE_REGION){
//		    qDebug() << iThetaI << iThetaS << iPhiS << smatrix[m_scatAnglesIndexSameDir[iThetaI][iThetaS][iPhiS]][0][0] << smatrix[m_scatAnglesIndexSameDir[iThetaI][iThetaS][iPhiS]][0][1] << smatrix[m_scatAnglesIndexSameDir[iThetaI][iThetaS][iPhiS]][1][1] << smatrix[m_scatAnglesIndexSameDir[iThetaI][iThetaS][iPhiS]][2][2] << "central value.";
		    double theta, phi, m11,  m12, m22, m33, sdaInt;
		    theta = phi = m11 = m12 = m22 = m33 = sdaInt = 0;
		    for(theta = thetaS - ANGLE_STEP/2; theta < thetaS + ANGLE_STEP/2; theta += STEP_AVG){
			for(phi = phiS-ANGLE_STEP/2; phi < phiS + ANGLE_STEP/2; phi += STEP_AVG){
			    double muScat, scatAngle;
			    muScat = Cos(thetaI)*Cos(theta+STEP_AVG/2) + Cos(phiI - (phi+STEP_AVG/2))*Sin(thetaI)*Sin(theta+STEP_AVG/2);
			    scatAngle = acos(muScat);
			    int index;
			    index = qRound(scatAngle / STEP_SMALL_SCAT);

			    m11 += smatrix[index][0][0] * Sin(theta+STEP_AVG/2) * STEP_AVG * STEP_AVG;
			    m12 += smatrix[index][0][1] * Sin(theta+STEP_AVG/2) * STEP_AVG * STEP_AVG;
			    m22 += smatrix[index][1][1] * Sin(theta+STEP_AVG/2) * STEP_AVG * STEP_AVG;
			    m33 += smatrix[index][2][2] * Sin(theta+STEP_AVG/2) * STEP_AVG * STEP_AVG;
			    sdaInt += Sin(theta+STEP_AVG/2) * STEP_AVG * STEP_AVG;
//			    if(iThetaI==1 && iThetaS==2 && iPhiS==48)
//				qDebug() << scatAngle << index << smatrix[index][0][0];
			}
		    }

		    pbsdfAero[iThetaI][iThetaS][iPhiS][0][0] = m11/sdaInt;
		    pbsdfAero[iThetaI][iThetaS][iPhiS][0][1] = m12/sdaInt;
		    pbsdfAero[iThetaI][iThetaS][iPhiS][1][1] = m22/sdaInt;
		    pbsdfAero[iThetaI][iThetaS][iPhiS][2][2] = m33/sdaInt;
//		    qDebug() << "averaged" << pbsdfAero[iThetaI][iThetaS][iPhiS][0][0] << pbsdfAero[iThetaI][iThetaS][iPhiS][0][1] << pbsdfAero[iThetaI][iThetaS][iPhiS][1][1] << pbsdfAero[iThetaI][iThetaS][iPhiS][2][2];
		}

	    }
	}
    }

}
